﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Runtime.Serialization;
using System.ServiceModel;
using System.Text;
using RaikesSimplexService.Contracts;
using RaikesSimplexService.DataModel;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.LinearAlgebra.Generic;
/*
 * Raikes Simplex Solver
 * SlyFox - Caitlin Bales, Derek Christensen, Travis Schreier, and Taylor Smith
 * This program solves simplex equations.
 */

namespace RaikesSimplexService.SlyFox
{

    public class Solver : ISolver
    {
        /// <summary>
        /// The Solve method solves the given simplex model
        /// </summary>
        /// <param name="model">The simplex problem to be solved</param>
        /// <returns>The solution to the simplex problem</returns>
        public Solution Solve(Model model)
        {
            List<LinearConstraint> constraints = model.Constraints;
            Solution solution = new Solution();
            Goal modelGoal = model.Goal;
            var rowVectors = new List<Vector<double>>(); 
            var costRow = model.Goal.Coefficients;

            //If the goal is a minimization problem, changes it to a maximization one.
            if (model.GoalKind == DataModel.GoalKind.Minimize)
            {
                for (int i = 0; i < model.Goal.Coefficients.Length; i++)
                    model.Goal.Coefficients[i] *= -1;
            }

            //Makes an array of coefficients, adding slack, surplus, or artificial variables as necessary
            var surplusCounter = 0;
            var artificialCounter = 0;
            var artificialList = new List<int>();
            var bv = new int [constraints.Count];
            List<int> basicVariables;

            for (int i = 0; i < constraints.Count; i++)
            {
                var constraint = (constraints.ElementAt(i).Coefficients);
                DenseVector constraintVector;
                 
                //Fills in the constraint vector with its current values and an artificial variable.
                if (constraints.ElementAt(i).Relationship == DataModel.Relationship.Equals)
                {
                    constraintVector = new DenseVector(constraint.Length+surplusCounter+i+1);
                    for (int j = 0; j < constraint.Length; j++)
                        constraintVector[j] = constraint[j];

                    constraintVector[surplusCounter+i+constraint.Length] = 1;
                    artificialList.Add(constraintVector.Count - 1);
                    bv[i]=constraintVector.Count -1;
                    artificialCounter++;
                    rowVectors.Add(constraintVector);
                }

                //Fills in the constraint vector with its current values and a slack variable.
                else if (constraints.ElementAt(i).Relationship == DataModel.Relationship.LessThanOrEquals)
                {
                    constraintVector = new DenseVector(constraint.Length + surplusCounter + i + 1);
                    for (int j = 0; j < constraint.Length; j++)
                        constraintVector[j] = constraint[j];

                    constraintVector[surplusCounter + i + constraint.Length] = 1;
                    bv[i]=constraintVector.Count-1;
                    rowVectors.Add(constraintVector);
                }

                //Fills in the constraint vector with its current values, a surplus variable, and an artificial variable.
                else if (constraints.ElementAt(i).Relationship == DataModel.Relationship.GreaterThanOrEquals)
                {
                    constraintVector = new DenseVector(constraint.Length + surplusCounter + i + 2);
                    for (int j = 0; j < constraint.Length; j++)
                        constraintVector[j] = constraint[j];

                    constraintVector[surplusCounter + i + constraint.Length] = -1;
                    constraintVector[surplusCounter + i + constraint.Length + 1] = 1;
                    surplusCounter++;
                    artificialList.Add(constraintVector.Count - 1);
                    artificialCounter++;
                    bv[i]=constraintVector.Count-1;
                    rowVectors.Add(constraintVector);
                }
            }

            Vector<double> RHSList = new DenseVector(constraints.Count);

            for (int i = 0; i < constraints.Count; i++)
            {
                RHSList[i] = model.Constraints.ElementAt(i).Value;
            }

            var b = RHSList.ToColumnMatrix();
            var denseMatrix = DenseMatrix.CreateFromRows(rowVectors);
            Vector<double> cRows = new DenseVector(constraints.Count + surplusCounter + costRow.Length);
            var CStar = new double[cRows.Count];

            //Runs only if we need to do phase one (there are artificial variables)
            if (artificialCounter > 0)
            {
                Vector<double> cRowsPhase1 = new DenseVector(constraints.Count + surplusCounter + costRow.Length);
                
                //make cRowsPhase1 (Min W)
                for (int i = 0; i < artificialList.Count; i++)
                {
                    cRowsPhase1[artificialList[i]] = -1;
                }

                //Runs phase one
                phase(ref denseMatrix, ref bv, ref cRowsPhase1, ref CStar, ref solution, ref b);

                basicVariables = bv.ToList();
                for (int i = 0; i < artificialList.Count; i++)
                {
                    //If we still have artificial basic variables, the solution is infeasible.
                    if(basicVariables.Contains(artificialList[i])){
                        solution.IsInfeasible = true;
                        return solution;
                    }
                }
            }

            cRows = new DenseVector(denseMatrix.ColumnCount - artificialCounter); 
            var columnsOfMatrixWithoutArtificialVariables = new List<Vector<double>>(denseMatrix.ColumnCount-artificialCounter);
            CStar = new double[cRows.Count];
            
            //removes artificial variables columns
            for (int i = 0; i < denseMatrix.ColumnCount; i++)
            {
                if(!artificialList.Contains(i))
                    columnsOfMatrixWithoutArtificialVariables.Add(denseMatrix.Column(i));
            }

           //Removes the artificial variables
           basicVariables = bv.ToList();
           int secondToLast = -1, last, decrement;

           for (int i = 0; i < denseMatrix.ColumnCount; i++)
           {
                if (!artificialList.Contains(i))
                {
                    last = i;
                    decrement = last - secondToLast - 1;
                    secondToLast = last - decrement;
                    if (basicVariables.Contains(i))
                        basicVariables[basicVariables.IndexOf(i)] = secondToLast;
                }
            }

            bv = basicVariables.ToArray();
            denseMatrix = DenseMatrix.CreateFromColumns(columnsOfMatrixWithoutArtificialVariables);

            //Updates the cost rows
            for (int i=0; i<costRow.Length; i++)
                cRows[i]=costRow[i];

            //Runs phase two on the matrix
            var bStar = phase(ref denseMatrix, ref bv, ref cRows, ref CStar, ref solution, ref b);

            basicVariables = bv.ToList();
            for (int i = 0; i < CStar.Length; i++)
            {
                if (equalToZero(CStar[i]) && !basicVariables.Contains(i))
                {
                    solution.AlternateSolutionsExist = true;
                }
            }

            var endValues = new double[costRow.Length];
            
            for (int i = 0; i < costRow.Length; i++)
            {
                if (bv.Contains(i))
                {
                    endValues[i] = bStar[basicVariables.IndexOf(i), 0];
                    endValues[i] = Math.Round(endValues[i], 3);
                }
            }
            if (model.GoalKind == DataModel.GoalKind.Minimize)
            {
                for (int i = 0; i < model.Goal.Coefficients.Length; i++)
                    model.Goal.Coefficients[i] *= -1;
            }

            solution.Values = endValues;
            return solution;
        }

        /// <summary>
        /// This creates the B-matrix, the matrix which contains
        /// the problem's basic variables.
        /// </summary>
        /// <param name="denseMatrix">The matrix of all constraints</param>
        /// <param name="basicVariables">The array of where the basic variables are held</param>
        /// <returns>A B-matrix</returns>
        public Matrix<double> createBMatrix(Matrix<double> denseMatrix, int[] basicVariables)
        {
            var bvPs = new List<Vector<double>>();
            for (int i = 0; i < basicVariables.Length; i++)
            {
                DenseVector P = new DenseVector(denseMatrix.Column(basicVariables[i]));
                bvPs.Add(P);
            }

            var BMatrix = DenseMatrix.CreateFromColumns(bvPs);
            return BMatrix;
        }

        /// <summary>
        /// Calculates P*
        /// </summary>
        /// <param name="inverseMatrix">The inverse of the B-matrix</param>
        /// <param name="originalP">The original P values</param>
        /// <returns>A P* vector</returns>
        public Matrix<double> findPStar(Matrix<double> inverseMatrix, Vector<double> originalP)
        {
            var P = originalP.ToColumnMatrix();
            var PStar = inverseMatrix.Multiply(P);
            return PStar;
        }

        /// <summary>
        /// Finds C*, the updated cost
        /// </summary>
        /// <param name="PStar">The updated constraint coefficients</param>
        /// <param name="basicVarOriginalC">The original cost coefficients</param>
        /// <param name="originalC">The original cost</param>
        /// <returns>C*, the updated cost</returns>
        public double findCStar(Matrix<double> PStar, Vector<double> basicVarOriginalC, double originalC)
        {
            var bvoc = basicVarOriginalC.ToRowMatrix();
            var product = bvoc.Multiply(PStar);
            double CStar = originalC - product[0,0];
            return CStar;
        }

        /// <summary>
        /// Finds b*, the updated right-hand side variables
        /// </summary>
        /// <param name="inverseMatrix">The inverted B-matrix</param>
        /// <param name="b">The original b matrix</param>
        /// <returns>The b* matrix</returns>
        public Matrix<double> findbStar(Matrix<double> inverseMatrix, Matrix<double> b)
        {
            var bStar = inverseMatrix.Multiply(b);
            return bStar;
        }

        /// <summary>
        /// Creates a tolerance for zero
        /// </summary>
        /// <param name="zero">The double to compare to zero</param>
        /// <returns>Whether the value is close enough to zero</returns>
        public Boolean equalToZero(double zero)
        {
            return (zero > -.00001 && zero < .00001);
        }

        /// <summary>
        /// Creates a tolerance for zero
        /// </summary>
        /// <param name="zero">The double to compare to zero</param>
        /// <returns>Whether the value is greater than zero</returns>
        public Boolean greaterThanZero(double zero)
        {
            return (zero > .00001);
        }

        /// <summary>
        /// Creates a tolerance for zero
        /// </summary>
        /// <param name="zero">The double to compare to zero</param>
        /// <returns>Whether the value is less than zero</returns>
        public Boolean lessThanZero(double zero)
        {
            return (zero < -.00001);
        }

        /// <summary>
        /// Runs a phase (either first or second) on the data.
        /// This method calculates the necessary matrices and pivots
        /// when it finds an updated cost that is greater than zero.
        /// It also checks to make sure that the solution is neither unbounded
        /// nor cycles.
        /// </summary>
        /// <param name="denseMatrix">The matrix of all constraints</param>
        /// <param name="bv">The array of basic variables</param>
        /// <param name="cRows">The vector of cost coefficients</param>
        /// <param name="CStar">The array of updated cost coefficients</param>
        /// <param name="solution">The problem's solution</param>
        /// <param name="b">The b matrix</param>
        /// <returns>The b* matrix</returns>
        public Matrix<double> phase(ref Matrix<double> denseMatrix, ref int[] bv, ref Vector<double> cRows, ref double[] CStar, ref Solution solution, ref Matrix<double> b)
        {
            Boolean timeToPivot = true;
            int timeoutCounter = 0;
            Boolean unbounded = false;
            var BMatrix = createBMatrix(denseMatrix, bv);
            var BInverse = BMatrix.Inverse();
            var bStar = findbStar(BInverse, b);
            //This will run until an optimal solution or exception is found.
            while (timeToPivot == true && timeoutCounter < 10000 && unbounded == false && solution.IsInfeasible==false)
            {
                BMatrix = createBMatrix(denseMatrix, bv);
                BInverse = BMatrix.Inverse();
                timeToPivot = false;
                int k = 0, leavingBV = -1;

                //Runs until a suitable pivot value is found.
                while (timeToPivot == false && k < denseMatrix.ColumnCount)
                {
                    if (!bv.Contains(k))
                    {
                        var cRowBV = new DenseVector(bv.Length);

                        for (int i = 0; i < bv.Length; i++)
                            cRowBV[i] = cRows[bv[i]];

                        var PRow = denseMatrix.Column(k);
                        var PStar = findPStar(BInverse, PRow);
                        CStar[k] = findCStar(PStar, cRowBV, cRows.ElementAt(k));

                        //If C* is greater than zero, it is a suitable pivot
                        if (greaterThanZero(CStar[k]))
                        {
                            timeToPivot = true;
                            bStar = findbStar(BInverse, b);
                            for (int i = 0; i < bStar.RowCount; i++)
                            {
                                if (bStar[i, 0] < 0)
                                    solution.IsInfeasible=true;
                            }
                            double min = -1;
                            Boolean minNotFound = true;
                            leavingBV = -1;

                            //Finds the minimum ratio of b* to P*
                            for (int i = 0; i < PRow.Count; i++)
                            {
                                double ratio = bStar[i, 0] / PStar[i,0];

                                if ((greaterThanZero(ratio) || equalToZero(ratio)) && minNotFound == true && PStar[i,0] > 0)
                                {
                                    min = ratio;
                                    leavingBV = i;
                                    minNotFound = false;
                                }

                                if ((greaterThanZero(ratio) || equalToZero(ratio)) && ratio < min && PStar[i,0] > 0) //SAME AS LAST IF STATEMENT???
                                {
                                    leavingBV = i;
                                    min = ratio;
                                }
                            }
                            //The solution is unbounded
                            if (min == -1)
                                unbounded = true;
                        }
                    }
                    else
                        CStar[k] = 0;

                    if (leavingBV != -1)
                        bv[leavingBV] = k;
                    k++;
                }
                timeoutCounter++;
               
            }
            //Allows a certain number of iterations before the calculation times out. Also
            //returns a timeout exception if it is unbounded.
            if (timeoutCounter >= 10000 || unbounded == true)
                solution.TimedOut = true;

            bStar = findbStar(BInverse, b);
            return bStar;
        }
    }
}

